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The Feynman integral is one of the most accurate methods for calculating density operator dy- 
namics in open quantum systems, but due to its enormous computational cost, one can only use it 
to calculate an approximation of the density operator at a sparse grid of points in time. Conven- 
tionally, interpolation methods such as splines are then used to to estimate the density operator in 
between these points, but this can lead to serious problems such as the loss of posit ivity. In this 
work a method is presented which uses physical information about the system to improve this inter- 
polation. The method is tested on a physically significant system and allows for a huge reduction 
in the amount of memory and CPU time required. 
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INTRODUCTION 



The statistical properties of a quantum state can be 
uniquely described by its density operator p. We are 
often interested in how the density operator of an open 
quantum system (OQS) evolves in time, and an enormous 
amount of effort has been dedicated towards developing 
mathematical techniques for predicting this evolution in 
time. 

To date, one of the most powerful of these techniques 
is based on the Feynman integral [2] and the Feynman- 
Vernon influence functional [2j [7]. This method, whose 
computation essentially only relies on integration, has 
substantial advantages over methods which rely on differ- 
ential equations for the density operator (quantum mas- 
ter equations). 

Unfortunately, although the Feynman integral meth- 
ods are conceptually straightforward, the cost of their 
numerical evaluation often increases exponentially in the 
number of time steps used. Contrarily, very small time 
steps are required in situations where p varies quickly. 

In fact, with current technology, using these methods 
is often intractable for the results desired. To reduce the 
computation time, one can store numbers that are reused 
in more than one iteration. But merely storing the am- 
plitude over each Feynman path ordinarily requires an 
array of M 2 ^ +1 ^ complex-valued elements (M = dimen- 
sion of the hilbert space of the OQS, or of the number of 
points in the discretization of the system coordinate in 
the case of an infinite dimensional hilbert space, and TV 
= number of time steps, excluding t = 0, elapsed before 
the influence functional becomes constant with respect 
to time). For a two-level system (2LS) with TV = 20, this 



amounts to > 70TE^]of memory (assuming double pre- 
cision arithmetic). If TV is increased to 25, this becomes 
> 72PE0 

The amount of memory required can be significantly 
reduced without too big of an impact on the accu- 
racy, by filtering out some of the smaller Feynman 
amplitudes [5j |6], but despite such filtering, there is al- 
ways a maximum TV (henceforth denoted TV max ) beyond 
which there is no longer enough memory for the surviving 
amplitudes to be stored, and the remaining ones would 
have to be recaluclated at each iteration (an excruciat- 
ingly time-demanding process). 

The point at which TV can no longer practically be in- 
creased, and still does not meet the demands of the situ- 
ation, occurs very often in real systems of interest. Even 
in the simple case of a driven 2LS, when the driving fre- 
quency becomes large enough, using the maximum com- 
putationally feasible TV will lead to an approximation of 
p at a grid of points which is too sparse for typical inter- 
polation techniques to capture important features such 
as peaks and troughs. 

This paper explores a method that attempts to capture 
these features, without increasing the amount of memory 
required. Instead of interpolating between these sparse 
points using conventional techniques, the interpolation 
incorporates some of the information about the relevant 
physics of the system. 



1 Throughout this paper, the SI convention will be used, so 
1TB=10 12 bytes rather than the binary convention in which 
lTB=lTiB=2 40 bytes. 

2 At the top of the current TOP500 list of supercomputers is the 
Cray XT5 Jaguar which has hard drive space on the order of 
lOPBpQ 
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EXAMPLE CALCULATIONS 

All example calculations below will be for a 2LS with 
hamiltonian: H = ^(|0)(1| + |1)(0|), and spectral den- 
sity: J(uj) = 0.0268 ^3 e (-o.2063a; 2 ) at a temperature of 
25K. This hamiltonian corresponds to a GaAs quantum 
dot being driven by a pulse with frequency ft = , and 
the spectral density is from a very recent experimental 
study [5]. The density matrix in every case will be ini- 
tialized at p = |0)(0|, and its evolution will be calculated 
using the QUAPI as described in [4] for a duration of 
3.5ps. For this entire duration, the influence functional 
for this spectral density never settles with respect to time, 
so the markovian iteration method described in [3] will 
not be used. 

I will also assume that the amplitude over each Feyn- 
man path is stored in an array named A, rather than 
being recalculated at each time step. For pedagogical 
reasons, no filtering of these amplitudes is performed in 
this paper, but the methods described here can easily be 
combined with a filtering process such as [6] or [5] to save 
computational cost. 



For such cases, no conventional interpolation technique 
would be able to capture qualitative features such as in 
the figure above. We need a better way to approximate 
p at places between the points of the Feynman mesh. 

To do this, let's construct a finer mesh, for example: 

i i | i i | I i | i i | i i 

Here, the blue (coarser) mesh is the Feynman mesh - 
for each point on this mesh, the size of A increases by a 
factor of M 2 . The ambition is to approximate p at the 
points on the red (finer) mesh, in a manner that not only 
resembles the true physics, but also does not increase the 
number of amplitudes over which to sum. 

There is a very simple way to do this. On the points of 
the Feynman mesh, p is still calculated exactly as before, 
so the value of p at these points is in no way damaged 
and this calculation can be completed before anything 
else. Then we can use a master equation to approximate 
p in between these points. A simple example is the von 
Neumann equation for a closed quantum system with the 
same system hamiltonian. If the hamiltonian is time- 
independent, this is: 



THE METHOD 

Conventionally, the numerical calculation of a Feyn- 
man integral involves dividing the time axis into a dis- 
crete mesh. As explained above, the number of points 
on this mesh is < some A^ max (this excludes the starting 
point t = 0). Let's define this mesh the Feynman mesh. 

Two calculations are shown below for the example sys- 
tem described above, each with a different number of 
time steps N. The square-shaped markers indicate the 
points on the Feynman mesh for each specific curve, and 
the points are joined by cubic splines (by the SPLINE 
function in MATLAB2010a). 
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Time (picoseconds) 



It's clear that with N — 2, some important qualita- 
tive features are missing. Unfortunately, some physically 
interesting quantum mechanical systems have a large 
enough hilbert space that using TV as low as 3 would 
currently be unachievable on a common computer. One 
such example is the hilbert space of single excitations in 
the the light-harvesting system LH2, the dimension of 
which is 27 [? ]. 



p(t + nAt) = U n p(t)U^ n , U = e~ iHA \ (1) 

where At is the size of the time-steps on the finer of 
the two meshes, and t is a point in time that is on the 
Feynman mesh. For time-dependent hamiltonians, this 
becomes: 

p(t + nAt) = e" 1 H ^ At 'p{ty V" *(*'><«' . ( 2 ) 

The figure below shows that when p is approximated 
at the points of the finer mesh in this way, the qualita- 
tive features of a QUAPI calculation with N = 12 can be 
recovered quite impressively by a calculation with only 
N = 2, provided that the number (NN) of points on the 
finer mesh lying between points of the Feynman mesh is 
increased to 2. The meshes for the upper diagram were 
designed such that the three points on the N2NN0 Feyn- 
man mesh also happen to be points on the Feynman mesh 
of all other cases (the total number of points on the mesh 
is 1 + (NN + 1)N). The meshes for the lower diagram 
use the maximum number of evenly spaced points that 
can fit in the given time interval (1 + (NN + 1)(N + 1)), 
so the size of the time steps for a given pair (TV, NN) is 
reduced in order to fit an extra NN points. The solid 
black curve was calculated with N12NN0, and is the ac- 
curate benchmark to which all other calculations in this 
paper are compared. The dotted black curve is what the 
evolution would look like if the OQS was closed (ie, uni- 
tary dynamics according to the von Neumann equation 
without any system-environment interaction). 
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In the upper diagram we can see that all curves with 
N = 2 have two points in common, which are the points 
on the Feynman mesh. This illustrates the fact (that was 
stated earlier), that the calculated value of p at the points 
of the Feynman mesh are in no way affected by this in- 
terpolation technique. The curves are simply clamped by 
the points on the Feynman mesh, and in between those 
points they evolve as if there was no interaction between 
the system and the environment. Because p is not ap- 
proximated very well at the points of the N = 2 Feynman 
mesh, all of the curves are clamped at somewhat inaccu- 
rate values, and therefore, regardless of how accurate the 
unitary approximation is on the finer mesh, some fea- 
tures such as the maxima and minima are not captured 
very well. 

In the lower diagram, the curves are clamped at differ- 
ent points because the size of the time steps on the Feyn- 
man mesh is reduced as NN increases. As a result, the 
maxima and minima are approximated slightly better. In 
fact, if we allowed the points on the Feynman mesh to be 
distributed unevenly, the curves could be clamped where 
the minima are predicted to be, and with N still equal to 
2, we could get results that look even better than those 
in the lower diagram. But to keep this study simple, all 
results in this paper will use meshes with evenly spaced 
points. 



FASTER CONVERGENCE 

When a Feynman integral is numerically approxi- 
mated, its accuracy is usually checked by systematically 



increasing the number of time steps and monitoring the 
results until they have converged. Once increasing the 
number of time steps leads to no change in the result (up 
to some tolerance), it is believed that the converged re- 
sult is accurate (provided that the error in all approxima- 
tions, such as the Trotter splitting of the time evolution 
operator, are only dependent on the size of the time step 
- this is usually true). 

The diagram below shows that for larger values of N A/", 
the curves converge quicker as N is increased. The size of 
the double-precision array A and the approximate CPU 
time for each calculation are shown in the legend - neither 
of these quantities changed at all as NN was varied. 
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ACCURACY 

We have looked at diagrams in which N was fixed and 
NN varied, and diagrams in which NN was fixed and N 
was varied. It is certainly apparent from fig. 2 that in- 
creasing NN allows calculations with low N to capture 
qualitative features that wouldn't be captured conven- 
tionally (with NN — 0). But since N is only 2, and 
therefore p is not clamped very accurately at the points 
on the Feynman mesh, this diagram doesn't convince us 
that a calculation with a high N can be accurately re- 
produced by using a lower N with a higher TV N. 

It is also certainly apparent from fig. 3 that by in- 
creasing NN, the calculations converge up to a reason- 
able tolerance with a lower N, but since each diagram 
uses a different NN, it is not obvious whether or not the 
curves in every case are converging to the (most reliable) 
converged result of the NN = diagram. 

The figure below shows the most reliable converged re- 
sult (7V12AW0), and compares it to results for 7V10AW1, 
NSNN1 and N6NN1. For the latter three curves, the 
dotted lines represent their corresponding NN = re- 
sults; and the NN = 2 and 3 curves were not displayed 
because for these N values, those curves look essentially 

the same as the NN — 1 curves. 
1 r 
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This figure shows that with NN = 1, the N12N0 re- 
sult can quite accurately be reproduced with only N = 6, 
and even more accurately with N = 8. More impor- 
tantly, it shows that calculations done the conventional 
way (NN = 0) for the same two cases differ very largely 
from the desired result. This outcome can be extremely 
useful in cases where the demands of the calculation can 
not be met by the computational resources available (a 
very common situation for numerical Feynman integrals) . 

For example, if one desires to predict the decoherence 
rate of an entangled state of two qubits (a quantity that's 
tremendously useful when choosing material parameters 
for building a quantum computer), the dimension of the 
Hilbert space is 4, so just storing the array of double pre- 
cision amplitudes for 10 time steps would require 281TB 
of memory if no filtering is done or if filtering doesn't 
help. This is pragmatically impossible today, but if we 
assume that the result in the above figure doesn't com- 
pletely deteriorate when applied to the considered sys- 



tem, a remarkably accurate approximation can be calcu- 
lated with only 6 time steps. The equivalent array would 
only require 4.3GB, and would therefore easily fit in the 
RAM of a fairly good modern laptop computer. 

One more observation from fig. 4 is that although both 
N = 10 curves lie closer to the converged result than the 
N8NN1 result, there is not much of a difference between 
the NN = and NN = 1 cases, simply because the 
NN = case already seems to have enough points. It is 
important to realize that increasing TV TV will not always 
improve the accuracy of a calculation, especially when 
the size of the time intervals in the Feynman mesh are 
already very small in comparison to the time scale on 
which p changes. 



DISCUSSION 

One feature of the method described above that is 
subideal is that during the unitary evolution of p which 
occurs between the points of the Feynman mesh, no in- 
formation about the influence of the environment on the 
OQS is used. This could very easily be improved - in- 
stead of using the solution of the von Neumann equation 
for a closed system to propagate p in these time inter- 
vals, one could use a more sophisticated master equation 
that incorporates information about the influence of the 
environment. 

Another feature that is independently subideal, is that 
even though p is being approximated at points in between 
the points of the Feynman mesh, no information is being 
used to try to update the amplitudes during these time 
intervals. Consequently, every approximation of p on a 
point of the Feynman mesh is only as accurate as it would 
have been (with a Feynman mesh of the same size) if this 
method was not used at all. 

For example, the interpolation technique described 
above neither improves nor impairs the accuracy of plots 
where the only points considered are on the Feynman 
mesh, such as ones where the elements of p(t = 15ps) are 
plotted as a function of driving frequency. If t = 15ps 
happens to be on the Feynman mesh, absolutely no im- 
provement is made in the approximation of p(t = 15ps). 
Even if the size of the Feynman time intervals were short- 
ened without changing N (in the same way they are 
shortened from the fig. 2a to fig. 2b), the improvement 
this method would provide would be negligible if t is large 
enough, which would usually be the case. 

There is more than one way in which physical infor- 
mation can be used to update the amplitudes more often 
without increasing the size of the array containing them. 
For one, we could continue to update p in a manner re- 
sembling eq. ([I]), but simultaneously use the operators 
in eq. ([!]) to update the elements of A. Alternatively, 
instead of using to update p on the points of the 
finer mesh, we could update these points by calculating 
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a Feynman integral, except instead of increasing the num- 
ber of amplitudes by a factor of M 2 (as would be done at 
each point on the Feynman mesh), the summation could 
be done on updated versions of the amplitudes left over 
from the previous point on the Feynman mesh. These up- 
dated amplitudes would ideally contain information from 
both the closed system hamiltonian and from the influ- 
ence functional. One could either use the influence func- 
tional from the previous Feynman point to do this, or 
one could actually update the influence functional along 
the way. 

These enhancements to the method described in this 
paper will be investigated much more thoroughly in fu- 
ture work. 
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